Regularized Partial Correlation Network
Adapted from Epskamp & Fried (2017)
Model selection example
Load packages
### 0. Load packages
if(!require(qgraph)) install.packages(pkgs="qgraph",repos="http://cran.r-project.org")## Loading required package: qgraph
## Registered S3 methods overwritten by 'huge':
## method from
## plot.sim BDgraph
## print.sim BDgraph
## Loading required package: bootnet
## Loading required package: ggplot2
## This is bootnet 1.4.3
## For questions and issues, please see github.com/SachaEpskamp/bootnet.
## Loading required package: dplyr
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Loading required package: igraph
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:dplyr':
##
## as_data_frame, groups, union
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
## Loading required package: doParallel
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel
Generate true network
### 1. Model selection example (figure 1 & 2)
# Set random seed:
set.seed(3)
# Generate true network:
# relationship between nodes
# an adjacency matrix: a square matrix used to represent a finite graph
Kappa <- as.matrix(get.adjacency(watts.strogatz.game(1,8,1,0)))
Kappa ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0 1 0 0 0 0 0 1
## [2,] 1 0 1 0 0 0 0 0
## [3,] 0 1 0 1 0 0 0 0
## [4,] 0 0 1 0 1 0 0 0
## [5,] 0 0 0 1 0 1 0 0
## [6,] 0 0 0 0 1 0 1 0
## [7,] 0 0 0 0 0 1 0 1
## [8,] 1 0 0 0 0 0 1 0
help("watts.strogatz.game")
The Watts-Strogatz small-world model
Description
Generate a graph according to the Watts-Strogatz network model.
Usage
sample_smallworld(dim, size, nei, p, loops = FALSE, multiple = FALSE)
smallworld(...)
Arguments
dim
Integer constant, the dimension of the starting lattice.
size
Integer constant, the size of the lattice along each dimension.
nei
Integer constant, the neighborhood within which the vertices of the lattice will be connected.
p
Real constant between zero and one, the rewiring probability.
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0 1 0 0 0 0 0 1
## [2,] 1 0 1 0 0 0 0 0
## [3,] 0 1 0 1 0 0 0 0
## [4,] 0 0 1 0 1 0 0 0
## [5,] 0 0 0 1 0 1 0 0
## [6,] 0 0 0 0 1 0 1 0
## [7,] 0 0 0 0 0 1 0 1
## [8,] 1 0 0 0 0 0 1 0
# coefficients between nodes (upper.tri)
Kappa[upper.tri(Kappa)] <- runif(sum(upper.tri(Kappa)),0.3,0.4) * Kappa[upper.tri(Kappa)]
# coefficients between nodes (upper.tri & lower.tri)
Kappa[lower.tri(Kappa)] <- t(Kappa)[lower.tri(Kappa)]
Kappa## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.0000000 0.3168042 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [2,] 0.3168042 0.0000000 0.3384942 0.0000000 0.0000000 0.0000000 0.0000000
## [3,] 0.0000000 0.3384942 0.0000000 0.3604394 0.0000000 0.0000000 0.0000000
## [4,] 0.0000000 0.0000000 0.3604394 0.0000000 0.3630979 0.0000000 0.0000000
## [5,] 0.0000000 0.0000000 0.0000000 0.3630979 0.0000000 0.3867919 0.0000000
## [6,] 0.0000000 0.0000000 0.0000000 0.0000000 0.3867919 0.0000000 0.3228202
## [7,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.3228202 0.0000000
## [8,] 0.3015330 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.3910148
## [,8]
## [1,] 0.3015330
## [2,] 0.0000000
## [3,] 0.0000000
## [4,] 0.0000000
## [5,] 0.0000000
## [6,] 0.0000000
## [7,] 0.3910148
## [8,] 0.0000000
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1.0000000 -0.3168042 0.0000000 0.0000000 0.0000000 0.0000000
## [2,] -0.3168042 1.0000000 -0.3384942 0.0000000 0.0000000 0.0000000
## [3,] 0.0000000 -0.3384942 1.0000000 -0.3604394 0.0000000 0.0000000
## [4,] 0.0000000 0.0000000 -0.3604394 1.0000000 -0.3630979 0.0000000
## [5,] 0.0000000 0.0000000 0.0000000 -0.3630979 1.0000000 -0.3867919
## [6,] 0.0000000 0.0000000 0.0000000 0.0000000 -0.3867919 1.0000000
## [7,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 -0.3228202
## [8,] -0.3015330 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [,7] [,8]
## [1,] 0.0000000 -0.3015330
## [2,] 0.0000000 0.0000000
## [3,] 0.0000000 0.0000000
## [4,] 0.0000000 0.0000000
## [5,] 0.0000000 0.0000000
## [6,] -0.3228202 0.0000000
## [7,] 1.0000000 -0.3910148
## [8,] -0.3910148 1.0000000
Model estimation
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1.00000000 0.5123636 0.14843729 -0.06473797 -0.03746480 -0.03790132
## [2,] 0.51236364 1.0000000 0.37934557 -0.23037365 -0.11652514 0.11409230
## [3,] 0.14843729 0.3793456 1.00000000 -0.42766411 -0.19930338 -0.06559148
## [4,] -0.06473797 -0.2303736 -0.42766411 1.00000000 0.52231864 0.21352558
## [5,] -0.03746480 -0.1165251 -0.19930338 0.52231864 1.00000000 0.33790299
## [6,] -0.03790132 0.1140923 -0.06559148 0.21352558 0.33790299 1.00000000
## [7,] 0.17674246 0.1262173 -0.06112605 0.19761716 0.20360528 0.36757685
## [8,] 0.48707457 0.2431745 0.03451387 0.15583973 0.06445881 0.06435186
## [,7] [,8]
## [1,] 0.17674246 0.48707457
## [2,] 0.12621730 0.24317454
## [3,] -0.06112605 0.03451387
## [4,] 0.19761716 0.15583973
## [5,] 0.20360528 0.06445881
## [6,] 0.36757685 0.06435186
## [7,] 1.00000000 0.39751551
## [8,] 0.39751551 1.00000000
# Labels:
Labs <- LETTERS[1:ncol(Kappa)]
### Figure 1: True network structure
# pdf("Fig1.pdf",width=4,height=4)
qgraph(qgraph:::wi2net(Kappa), labels = Labs, vsize = 15, edge.labels=TRUE, vTrans = 254, edge.label.cex=1.5, theme = "colorblind")# dev.off()
# EBIC analysis:
res0 <- EBICglasso(corMat,nrow(Data),0,returnAllResults=TRUE,nlambda=10)
res0.25 <- EBICglasso(corMat,nrow(Data),0.25,returnAllResults=TRUE,nlambda=10)
res0.5 <- EBICglasso(corMat,nrow(Data),0.5,returnAllResults=TRUE,nlambda=10)
### Figure 2: Estimated networks:
# pdf("Fig2.pdf",width=5*2,height=2*2)
layout(rbind(1:5,6:10))
for (i in 1:10){
qgraph(qgraph::wi2net(res0$results$wi[,,i]), labels = Labs, vsize = 15,mar=c(2,2,2,2), theme = "colorblind")
text(0,0.3,paste0("Lambda: ",round(res0$lambda[i],3)),cex=1)
text(0,0.1,paste0("EBIC (gamma = 0): ",round(res0$ebic[i],1)),cex=ifelse(res0$ebic[i] == min(res0$ebic),1,0.6),font = ifelse(res0$ebic[i] == min(res0$ebic),2,1))
text(0,-0.1,paste0("EBIC (gamma = 0.25): ",round(res0.25$ebic[i],1)),cex=ifelse(res0.25$ebic[i] == min(res0.25$ebic),1,0.6),font = ifelse(res0.25$ebic[i] == min(res0.25$ebic),2,1))
text(0,-0.3,paste0("EBIC (gamma = 0.5): ",round(res0.5$ebic[i],1)),cex=ifelse(res0.5$ebic[i] == min(res0.5$ebic),1,0.6),font = ifelse(res0.5$ebic[i] == min(res0.5$ebic),2,1))
}Conditional Element Selection
Description
ifelse returns a value with the same shape as test which is filled with elements selected from either yes or no depending on whether the element of test is TRUE or FALSE.
Usage
ifelse(test, yes, no)
Arguments
test
an object which can be coerced to logical mode.
yes
return values for true elements of test.
no
return values for false elements of test.
Empirical example
##### 2 DATA
# Read in spss data file (here from my own desktop, but you can easily change that of course)
data_full <- read.spss("https://www.lijinzhang.xyz/share/lab_pre/PTSDdata2.sav", to.data.frame=TRUE)## Warning in read.spss("https://www.lijinzhang.xyz/share/lab_pre/PTSDdata2.sav", :
## Undeclared level(s) 21, 25, 26, 27, 28, 29, 30, 31, 32, 35, 37, 40, 41, 43, 44,
## 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,
## 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 78, 81, 82, 84, 85, 86, 89 added
## in variable: PPAGE
# Create data frame of variables that we need for analysis:
## age, gender, all 20 symptoms (_MONTH extension), Sum_GAD2, Sum_PHQ2, Passive_SI_FINAL, Active_SI_FINAL, PCS, MCS, QualityofLife_SUM)
data <- as.data.frame(data_full[,c(11:30)]) # data for analysis
colnames(data) <- c(1:20)
head(data)## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 1 0 1 2 1 2 2 0 1 1 2 1 2 1 1 1 0 0 2 1
## 2 2 3 2 2 2 3 2 0 2 3 3 3 2 1 1 1 3 2 1 2
## 3 1 0 0 1 0 2 1 1 1 1 1 2 2 1 1 1 1 1 2 2
## 4 2 2 0 0 1 3 4 1 4 0 1 4 4 4 4 0 4 0 0 2
## 5 2 1 0 2 2 1 0 1 1 1 2 1 3 3 2 0 1 0 4 4
## 6 4 3 3 4 3 2 2 0 4 2 4 4 3 3 4 0 4 2 2 4
Model estimation & plot networks
##### 3 ESTIMATE & PLOT NETWORKS (Figure 3)
graph1 <- estimateNetwork(
data,
default = "EBICglasso",
corMethod = "cor_auto",
tuning = 0)## Estimating Network. Using package::function:
## - qgraph::EBICglasso for EBIC model selection
## - using glasso::glasso
## - qgraph::cor_auto for correlation computation
## - using lavaan::lavCor
## Variables detected as ordinal: 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20
## Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
## penalize.diagonal, : A dense regularized network was selected (lambda < 0.1 *
## lambda.max). Recent work indicates a possible drop in specificity. Interpret the
## presence of the smallest edges with care. Setting threshold = TRUE will enforce
## higher specificity, at the cost of sensitivity.
#default ="EBICglasso" specifies that the EBICglasso function from qgraph is used
graph2 <- estimateNetwork(
data,
default = "EBICglasso",
corMethod = "cor_auto",
tuning = 0.25)## Estimating Network. Using package::function:
## - qgraph::EBICglasso for EBIC model selection
## - using glasso::glasso
## - qgraph::cor_auto for correlation computation
## - using lavaan::lavCor
## Variables detected as ordinal: 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20
## Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
## penalize.diagonal, : A dense regularized network was selected (lambda < 0.1 *
## lambda.max). Recent work indicates a possible drop in specificity. Interpret the
## presence of the smallest edges with care. Setting threshold = TRUE will enforce
## higher specificity, at the cost of sensitivity.
## Estimating Network. Using package::function:
## - qgraph::EBICglasso for EBIC model selection
## - using glasso::glasso
## - qgraph::cor_auto for correlation computation
## - using lavaan::lavCor
## Variables detected as ordinal: 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20
#The corMethod = "cor_auto" argument specifies polychoric and polyserial correlations
# Compute average layout:
# To compare two networks, one should constrain the layout to be equal for both networks. One way to do so is by using averageLayout
L <- averageLayout(graph1, graph2, graph3)
#pdf("Fig3.pdf")
# maximum = 0.54 so it is consistent with later graphs
plot(graph1, layout=L, maximum = 0.54, cut = 0) ## 1 2 3 4 5 6
## 1 0.000000000 0.18874500 0.02041634 0.26949321 0.09830280 0.00000000
## 2 0.188745000 0.00000000 0.41772177 0.00000000 0.02984737 0.00000000
## 3 0.020416337 0.41772177 0.00000000 0.21831074 0.13159891 0.00000000
## 4 0.269493205 0.00000000 0.21831074 0.00000000 0.28643618 0.04458847
## 5 0.098302796 0.02984737 0.13159891 0.28643618 0.00000000 0.18003268
## 6 0.000000000 0.00000000 0.00000000 0.04458847 0.18003268 0.00000000
## 7 0.000000000 0.06858381 0.02724846 0.12789258 0.12053477 0.21212668
## 8 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.09115236
## 9 0.000000000 -0.12634394 0.00000000 0.03406929 0.00000000 -0.08216170
## 10 0.000000000 0.00000000 0.00000000 0.09573875 0.00654187 0.03321589
## 11 0.194457214 0.03799185 0.02511442 0.06800126 0.01263520 0.00000000
## 12 0.062786336 0.03125524 0.09709947 -0.03650199 0.00000000 0.00000000
## 13 0.000000000 0.00000000 0.00000000 -0.03771817 0.00000000 0.00000000
## 14 0.000000000 0.00000000 0.00000000 -0.10174362 0.00000000 0.00000000
## 15 0.002444043 0.00000000 0.00000000 0.00000000 0.11178363 0.00000000
## 16 -0.004295477 0.14300460 0.00000000 -0.04662989 0.00000000 0.00000000
## 17 0.000000000 0.00000000 0.07662094 0.00000000 0.00000000 0.06282331
## 18 -0.026095863 0.07743453 0.08133495 0.00000000 0.11272345 0.00000000
## 19 0.025083092 0.00000000 0.00000000 -0.01111239 0.00000000 0.00000000
## 20 0.005666384 0.11676674 0.09682865 0.00000000 0.03893225 0.00000000
## 7 8 9 10 11 12
## 1 0.00000000 0.00000000 0.00000000 0.00000000 0.19445721 0.06278634
## 2 0.06858381 0.00000000 -0.12634394 0.00000000 0.03799185 0.03125524
## 3 0.02724846 0.00000000 0.00000000 0.00000000 0.02511442 0.09709947
## 4 0.12789258 0.00000000 0.03406929 0.09573875 0.06800126 -0.03650199
## 5 0.12053477 0.00000000 0.00000000 0.00654187 0.01263520 0.00000000
## 6 0.21212668 0.09115236 -0.08216170 0.03321589 0.00000000 0.00000000
## 7 0.00000000 0.03972146 0.00000000 0.00000000 0.00000000 0.14892145
## 8 0.03972146 0.00000000 0.01061640 0.08314477 0.02193945 0.00000000
## 9 0.00000000 0.01061640 0.00000000 0.01903501 0.27003465 0.17438501
## 10 0.00000000 0.08314477 0.01903501 0.00000000 0.47855705 -0.03917531
## 11 0.00000000 0.02193945 0.27003465 0.47855705 0.00000000 0.00000000
## 12 0.14892145 0.00000000 0.17438501 -0.03917531 0.00000000 0.00000000
## 13 0.00000000 0.00000000 0.04846113 0.00000000 0.01398274 0.13173568
## 14 0.03207787 0.00000000 0.13011228 -0.02027354 0.00000000 0.09084672
## 15 0.08373357 0.01699078 0.10368434 0.00000000 0.16221888 0.06550974
## 16 0.06201111 0.11061140 0.00000000 0.05790148 0.00000000 0.00000000
## 17 0.00000000 0.05797455 0.05298531 0.00000000 0.00000000 0.00000000
## 18 0.00000000 0.00000000 0.00000000 -0.03159653 0.00000000 0.01513076
## 19 -0.11817257 0.10797522 0.05759475 0.01598574 0.02382759 0.30016698
## 20 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## 13 14 15 16 17 18
## 1 0.00000000 0.00000000 0.002444043 -0.004295477 0.00000000 -0.02609586
## 2 0.00000000 0.00000000 0.000000000 0.143004603 0.00000000 0.07743453
## 3 0.00000000 0.00000000 0.000000000 0.000000000 0.07662094 0.08133495
## 4 -0.03771817 -0.10174362 0.000000000 -0.046629892 0.00000000 0.00000000
## 5 0.00000000 0.00000000 0.111783627 0.000000000 0.00000000 0.11272345
## 6 0.00000000 0.00000000 0.000000000 0.000000000 0.06282331 0.00000000
## 7 0.00000000 0.03207787 0.083733569 0.062011109 0.00000000 0.00000000
## 8 0.00000000 0.00000000 0.016990779 0.110611398 0.05797455 0.00000000
## 9 0.04846113 0.13011228 0.103684338 0.000000000 0.05298531 0.00000000
## 10 0.00000000 -0.02027354 0.000000000 0.057901477 0.00000000 -0.03159653
## 11 0.01398274 0.00000000 0.162218881 0.000000000 0.00000000 0.00000000
## 12 0.13173568 0.09084672 0.065509741 0.000000000 0.00000000 0.01513076
## 13 0.00000000 0.47295317 0.025273362 0.061837398 0.08214261 0.09083379
## 14 0.47295317 0.00000000 0.113112038 0.000000000 0.00000000 0.00000000
## 15 0.02527336 0.11311204 0.000000000 0.266780501 0.00000000 0.00000000
## 16 0.06183740 0.00000000 0.266780501 0.000000000 0.03665648 0.11716494
## 17 0.08214261 0.00000000 0.000000000 0.036656477 0.00000000 0.32326412
## 18 0.09083379 0.00000000 0.000000000 0.117164939 0.32326412 0.00000000
## 19 0.15079596 0.02959822 0.012002099 0.000000000 0.00000000 0.09671241
## 20 0.11331185 0.02338779 0.000000000 0.000000000 0.01003874 0.10200155
## 19 20
## 1 0.02508309 0.005666384
## 2 0.00000000 0.116766744
## 3 0.00000000 0.096828655
## 4 -0.01111239 0.000000000
## 5 0.00000000 0.038932251
## 6 0.00000000 0.000000000
## 7 -0.11817257 0.000000000
## 8 0.10797522 0.000000000
## 9 0.05759475 0.000000000
## 10 0.01598574 0.000000000
## 11 0.02382759 0.000000000
## 12 0.30016698 0.000000000
## 13 0.15079596 0.113311854
## 14 0.02959822 0.023387794
## 15 0.01200210 0.000000000
## 16 0.00000000 0.000000000
## 17 0.00000000 0.010038744
## 18 0.09671241 0.102001553
## 19 0.00000000 0.201302918
## 20 0.20130292 0.000000000
## [1] 105
## [1] 95
## [1] 87
## [1] 5.091125
## [1] 4.625315
## [1] 4.337963
##### 4 CENTRALITY (Figure 4)
#pdf("Fig4.pdf")
centralityPlot(list(EBIC0=graph1,EBIC0.25=graph2,EBIC0.5=graph3), include="all")## Note: z-scores are shown on x-axis rather than raw centrality indices.
#dev.off()
# node strength, quantifying how well a node is directly connected to other nodes (该节点直接相连的连线加权值的绝对值之和)
# closeness, quantifying how well a node is indirectly connected to other nodes (所有节点到该节点的最短距离之和的倒数)
# betweenness, quantifying how important a node is in the average path between two other nodes. (该节点出现在网络中任意两个节点的最短路径上的次数)Prior sample size analysis
trueNetwork <- estimateNetwork(
data,
default = "EBICglasso",
corMethod = "cor_auto",
tuning = 0.5,
refit = TRUE)## Estimating Network. Using package::function:
## - qgraph::EBICglasso for EBIC model selection
## - using glasso::glasso
## - qgraph::cor_auto for correlation computation
## - using lavaan::lavCor
## Variables detected as ordinal: 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20
## Refitting network without LASSO regularization
library("bootnet")
simRes <- netSimulator(trueNetwork$graph,
dataGenerator = ggmGenerator(ordinal = TRUE, nLevels = 5),
default = "EBICglasso",
nCases = c(100,250,500,1000,2500),
tuning = 0.5,
nReps = 100,
nCores = 1
)## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
“# An empty network was selected to be the best fitting network. Possibly set ‘lambda.min.ratio’ higher to search more sparse networks. You can also change the ‘gamma’ parameter to improve sensitivity (at the cost of specificity).”
“## Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal = ## penalize.diagonal, : A dense regularized network was selected (lambda < 0.1 * ## lambda.max). Recent work indicates a possible drop in specificity. Interpret the ## presence of the smallest edges with care. Setting threshold = TRUE will enforce ## higher specificity, at the cost of sensitivity.”
Sensitivity: Also termed the true-positive rate, the proportion of edges present in the true network that were detected in the estimated network.
Specificity: Also termed the true-negative rate, the proportion of missing edges in the true network that were also detected correctly to be absent edges in the estimated network.
The correlation between edge weights of the true network and edge weights of the estimated network, or between centrality estimates based on the true network and centrality estimates based on the estimated network
## 1 2 3 4 5 6
## 1 0.000000000 0.195283249 0.005102409 0.27464323 0.093370806 0.00000000
## 2 0.195283249 0.000000000 0.457258512 0.00000000 0.017356395 0.00000000
## 3 0.005102409 0.457258512 0.000000000 0.21631094 0.141626695 0.00000000
## 4 0.274643230 0.000000000 0.216310940 0.00000000 0.303173116 0.03679962
## 5 0.093370806 0.017356395 0.141626695 0.30317312 0.000000000 0.18185716
## 6 0.000000000 0.000000000 0.000000000 0.03679962 0.181857157 0.00000000
## 7 0.000000000 0.067998717 0.015047013 0.12339396 0.122698337 0.22094378
## 8 0.000000000 0.000000000 0.000000000 0.00000000 0.000000000 0.09662244
## 9 0.000000000 0.000000000 0.000000000 0.00000000 0.000000000 0.00000000
## 10 0.000000000 0.000000000 0.000000000 0.09450993 -0.010580349 0.01414188
## 11 0.205659921 -0.010705344 0.021628378 0.07974091 0.007534505 0.00000000
## 12 0.064983356 -0.008753877 0.098501073 0.00000000 0.000000000 0.00000000
## 13 0.000000000 0.000000000 0.000000000 0.00000000 0.000000000 0.00000000
## 14 0.000000000 0.000000000 0.000000000 -0.20729516 0.000000000 0.00000000
## 15 -0.005972199 0.000000000 0.000000000 0.00000000 0.118998463 0.00000000
## 16 0.000000000 0.130339925 0.000000000 0.00000000 0.000000000 0.00000000
## 17 0.000000000 0.000000000 0.076840447 0.00000000 0.000000000 0.05562034
## 18 0.000000000 0.064252350 0.072374885 0.00000000 0.101640680 0.00000000
## 19 0.000000000 0.000000000 0.000000000 0.00000000 0.000000000 0.00000000
## 20 0.007479995 0.113931755 0.096802652 0.00000000 0.036211153 0.00000000
## 7 8 9 10 11 12
## 1 0.00000000 0.00000000 0.00000000 0.00000000 0.205659921 0.064983356
## 2 0.06799872 0.00000000 0.00000000 0.00000000 -0.010705344 -0.008753877
## 3 0.01504701 0.00000000 0.00000000 0.00000000 0.021628378 0.098501073
## 4 0.12339396 0.00000000 0.00000000 0.09450993 0.079740908 0.000000000
## 5 0.12269834 0.00000000 0.00000000 -0.01058035 0.007534505 0.000000000
## 6 0.22094378 0.09662244 0.00000000 0.01414188 0.000000000 0.000000000
## 7 0.00000000 0.03033601 0.00000000 0.00000000 0.000000000 0.106636372
## 8 0.03033601 0.00000000 0.00000000 0.09123952 0.022549418 0.000000000
## 9 0.00000000 0.00000000 0.00000000 0.00000000 0.271642439 0.146712454
## 10 0.00000000 0.09123952 0.00000000 0.00000000 0.505686959 0.000000000
## 11 0.00000000 0.02254942 0.27164244 0.50568696 0.000000000 0.000000000
## 12 0.10663637 0.00000000 0.14671245 0.00000000 0.000000000 0.000000000
## 13 0.00000000 0.00000000 0.03047257 0.00000000 0.000000000 0.129588198
## 14 0.00000000 0.00000000 0.15097476 0.00000000 0.000000000 0.118422046
## 15 0.08505491 0.02048741 0.07485573 0.00000000 0.180410481 0.067038344
## 16 0.04628658 0.12054616 0.00000000 0.03115945 0.000000000 0.000000000
## 17 0.00000000 0.06895300 0.02888292 0.00000000 0.000000000 0.000000000
## 18 0.00000000 0.00000000 0.00000000 0.00000000 0.000000000 0.010957557
## 19 0.00000000 0.10858793 0.04819525 0.00000000 0.024906667 0.299685751
## 20 0.00000000 0.00000000 0.00000000 0.00000000 0.000000000 0.000000000
## 13 14 15 16 17 18
## 1 0.00000000 0.00000000 -0.005972199 0.00000000 0.000000000 0.00000000
## 2 0.00000000 0.00000000 0.000000000 0.13033993 0.000000000 0.06425235
## 3 0.00000000 0.00000000 0.000000000 0.00000000 0.076840447 0.07237489
## 4 0.00000000 -0.20729516 0.000000000 0.00000000 0.000000000 0.00000000
## 5 0.00000000 0.00000000 0.118998463 0.00000000 0.000000000 0.10164068
## 6 0.00000000 0.00000000 0.000000000 0.00000000 0.055620338 0.00000000
## 7 0.00000000 0.00000000 0.085054913 0.04628658 0.000000000 0.00000000
## 8 0.00000000 0.00000000 0.020487413 0.12054616 0.068953005 0.00000000
## 9 0.03047257 0.15097476 0.074855727 0.00000000 0.028882922 0.00000000
## 10 0.00000000 0.00000000 0.000000000 0.03115945 0.000000000 0.00000000
## 11 0.00000000 0.00000000 0.180410481 0.00000000 0.000000000 0.00000000
## 12 0.12958820 0.11842205 0.067038344 0.00000000 0.000000000 0.01095756
## 13 0.00000000 0.50456045 0.011820020 0.06390397 0.089867535 0.09355390
## 14 0.50456045 0.00000000 0.153702903 0.00000000 0.000000000 0.00000000
## 15 0.01182002 0.15370290 0.000000000 0.27661091 -0.002980352 0.00000000
## 16 0.06390397 0.00000000 0.276610912 0.00000000 0.036502063 0.11124936
## 17 0.08986753 0.00000000 -0.002980352 0.03650206 0.000000000 0.34457546
## 18 0.09355390 0.00000000 0.000000000 0.11124936 0.344575459 0.00000000
## 19 0.15698509 0.02847355 -0.014459351 0.00000000 0.000000000 0.08224302
## 20 0.10422255 0.04081857 0.000000000 0.00000000 0.007031625 0.10910560
## 19 20
## 1 0.00000000 0.007479995
## 2 0.00000000 0.113931755
## 3 0.00000000 0.096802652
## 4 0.00000000 0.000000000
## 5 0.00000000 0.036211153
## 6 0.00000000 0.000000000
## 7 0.00000000 0.000000000
## 8 0.10858793 0.000000000
## 9 0.04819525 0.000000000
## 10 0.00000000 0.000000000
## 11 0.02490667 0.000000000
## 12 0.29968575 0.000000000
## 13 0.15698509 0.104222550
## 14 0.02847355 0.040818573
## 15 -0.01445935 0.000000000
## 16 0.00000000 0.000000000
## 17 0.00000000 0.007031625
## 18 0.08224302 0.109105596
## 19 0.00000000 0.201752964
## 20 0.20175296 0.000000000
## Warning: `select_()` is deprecated as of dplyr 0.7.0.
## Please use `select()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: `summarise_each_()` is deprecated as of dplyr 0.7.0.
## Please use `across()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## === netSimulator Results ===
##
## Mean (SD) values per varied levels:
##
## nCases default tuning ExpectedInfluence MaxFalseEdgeWidth bias
## 1 100 EBICglasso 0.5 0.8500182 0.3011283 0.04014474
## 2 500 EBICglasso 0.5 0.9213254 0.2011513 0.02145479
## 3 1000 EBICglasso 0.5 0.9549089 0.1258747 0.01717972
## 4 1000 EBICglasso 0.5 0.9714165 0.1874430 0.01843795
## 5 1000 EBICglasso 0.5 0.9763828 0.1581449 0.01749627
## 6 2500 EBICglasso 0.5 0.9818661 0.1106360 0.01429199
## sensitivity specificity correlation strength closeness betweenness
## 1 0.62 (NA) 0.72 (NA) 0.71 (NA) 0.81 (NA) 0.42 (NA) 0.49 (NA)
## 2 0.89 (NA) 0.62 (NA) 0.94 (NA) 0.93 (NA) 0.85 (NA) 0.92 (NA)
## 3 0.86 (NA) 0.68 (NA) 0.96 (NA) 0.95 (NA) 0.84 (NA) 0.87 (NA)
## 4 0.89 (NA) 0.64 (NA) 0.96 (NA) 0.98 (NA) 0.87 (NA) 0.86 (NA)
## 5 0.92 (NA) 0.65 (NA) 0.96 (NA) 0.97 (NA) 0.94 (NA) 0.93 (NA)
## 6 0.94 (NA) 0.51 (NA) 0.98 (NA) 0.99 (NA) 0.94 (NA) 0.8 (NA)
##
##
## Use plot(x) to plot results (nCases only), or as.data.frame(x) to see all results.
# N=250 achieves a correlation between the “true” and estimated networks above 0.8 for edge weights and strength, and above 0.7 for sensitivity.
plot(simRes,yvar = c("strength","closeness","betweenness"))"We recommend:
more research estimating network models from psychological data will make clear what one could expect as a true network structure, especially if researchers make the statistical parameters of their network models publicly available;
researchers should simulate network models under a wide variety of potential true network structures, using different estimation methods;
researchers should simulate data under an expected network structure to gain some insight in the required sample size. "
Post-hoc Stability checks
ncores <- 1
if(detectCores()-1 > 1) {
ncores <- detectCores()-5 # modify
}
boot1 <- try(bootnet(results, nCores = ncores, nBoots = 1000, type = "nonparametric"))## Note: bootnet will store only the following statistics: edge, strength, outStrength, inStrength
## Error in is(data, "bootnetResult") : 找不到对象'results'
How to debug?
## Note: bootnet will store only the following statistics: edge, strength, outStrength, inStrength
## Estimating sample network...
## Warning in formals(fun): argument is not a function
## Warning in formals(fun): argument is not a function
## Error in do.call(.input$estimator, c(list(data), .input$arguments)) :
## 'what' must be a function or character string
Search on the github: solution
graph3 <- estimateNetwork(
data,
default = "EBICglasso",
corMethod = "cor_auto",
tuning = 0.5)
boot1 <- bootnet(graph3, nCores = ncores, nBoots = 1000, type = "nonparametric")
boot2 <- bootnet(graph3, nCores = ncores, nBoots = 1000, type = "case", statistics=c("strength", "closeness","betweenness"))
# nonparametric bootstrap (using resampled data with replacement)
# Confidence intervals can not be constructed for centrality indices (see the supplementary materials of Epskamp, Borsboom et al., 2017).
# To assess the stability of centrality indices, one can perform a case-dropping bootstrap (subsampling without replacement).
## Bootstrap plots:
#pdf("BootnetResults.pdf", height = 10, width = 15)
plot(boot1, order = "sample")## Warning: `filter_()` is deprecated as of dplyr 0.7.0.
## Please use `filter()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: `summarise_()` is deprecated as of dplyr 0.7.0.
## Please use `summarise()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: `mutate_()` is deprecated as of dplyr 0.7.0.
## Please use `mutate()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
Figure S1. Bootstrapped edge-weights of the γ = 0.5 PTSD network shown in Figure 4.
## Expected significance level given number of bootstrap samples is approximately: 0.05
Figure S2. Bootstrapped significance (α = 0.05) between edges of the γ = 0.5 PTSD network shown in Figure 4. Eeach row and column indicates an edge. Black boxes represent significant differences and gray boxes represent non-significant differences. The color in the diagonal corresponds with the edge color shown in
## Expected significance level given number of bootstrap samples is approximately: 0.051
Figure S3. Bootstrapped significance (α = 0.05) between strength centrality metric of nodes of the γ = 0.5 PTSD network shown in Figure 4. Eeach row and column indicates a node. Black boxes represent significant differences and gray boxes represent non-significant differences. The value in the diagonal corresponds with the strength of a node.
Figure S4. The correlation between the original centrality index and the centrality index after dropping a percentage of subjects at random.
## === Correlation Stability Analysis ===
##
## Sampling levels tested:
## nPerson Drop% n
## 1 55 75.1 17
## 2 72 67.4 82
## 3 90 59.3 95
## 4 107 51.6 135
## 5 124 43.9 130
## 6 141 36.2 116
## 7 158 28.5 103
## 8 176 20.4 128
## 9 193 12.7 96
## 10 210 5.0 98
##
## Maximum drop proportions to retain correlation of 0.7 in at least 95% of the samples:
##
## betweenness: 0.05 (CS-coefficient is lowest level tested)
## - For more accuracy, run bootnet(..., caseMin = 0, caseMax = 0.127)
##
## closeness: 0.127
## - For more accuracy, run bootnet(..., caseMin = 0.05, caseMax = 0.204)
##
## strength: 0.516
## - For more accuracy, run bootnet(..., caseMin = 0.439, caseMax = 0.593)
##
## Accuracy can also be increased by increasing both 'nBoots' and 'caseN'.
"The results of the case-dropping bootstrap can also be summarized in a coefficient, the CS-coefficient (correlation stability), which quantifies the proportion of data that can be dropped to retain with 95% certainty a correlation of at least 0.7 with the original centrality coefficients.
Ideally this coefficient should be above 0.5, and should be at least above 0.25."
“Strength was shown to be stable (CS(cor=0.7) 0.516) while closeness (CS(cor=0.7) 0.204) and betweenness (CS(cor=0.7) 0.05) were not. Thus, the post hoc analysis shows that the estimated network structure and derived centrality indices should be interpreted with some care for our example network of PTSD symptoms.”